library(data.table)
library(lfe)
library(stargazer)
library(starpolishr)
library(ggplot2)
library(TwoWayFEWeights)
library(DIDmultiplegt)


rm(list=ls())
setwd("C:/Users/k2258581/OneDrive - King's College London/Drive/Research/Partisan Motivated Reasoning/Replication package")

############## Initialisation #############
iWindowLength <- 15
vMinYear_all <- 1974:2006 #starting year for shifting/shrinking window analyses

#import data
fData <- fread("GSS_allyears_clean.csv")



#import weights 2021 (not in original data)
fData_2021 <- readRDS("GSS_allyears_new.rds")
fData_2021 <- as.data.table(fData_2021)[year==2021][,c("year", "id","wtssps","realinc")]
fData <- fData[year==2021,wtssall := fData_2021$wtssps]


fData <- fData[!is.na(wtssall)] #remove data without weights
fData <- fData[year>=1974] #Correct data from 1974



#clean party identification
fData <- fData[,Party_cat := ifelse(partyid%in%c("independent, close to democrat",
                                                 "not very strong democrat",
                                                 "strong democrat"),"Democrat",
                                    ifelse(partyid%in%c("independent, close to republican",
                                                        "not very strong republican",
                                                        "strong republican"),"Republican",
                                           ifelse(partyid=="independent (neither, no response)","Independent",NA)))]

fData <- fData[,Dem := Party_cat=="Democrat"]
fData <- fData[,Rep := Party_cat=="Republican"]
fData <- fData[,Ind := Party_cat=="Independent"]

#define strong partisanship
fData <- fData[,Party_cat_strong := ifelse(partyid%in%c("strong democrat"),"Democrat",
                                           ifelse(partyid%in%c("strong republican"),"Republican",
                                                  ifelse(partyid=="independent (neither, no response)","Independent",NA)))]

fData <- fData[,Dem_strong := Party_cat_strong=="Democrat"]
fData <- fData[,Rep_strong := Party_cat_strong=="Republican"]
fData <- fData[,Ind_strong := Party_cat_strong=="Independent"]

#define semi-strong partisanship
fData <- fData[,Party_cat_semistrong := ifelse(partyid%in%c("strong democrat","not very strong democrat"),"Democrat",
                                               
                                               ifelse(partyid%in%c("strong republican","not very strong republican"),"Republican",
                                                      ifelse(partyid=="independent (neither, no response)","Independent",NA)))]

fData <- fData[,Dem_semistrong := Party_cat_semistrong=="Democrat"]
fData <- fData[,Rep_semistrong := Party_cat_semistrong=="Republican"]
fData <- fData[,Ind_semistrong := Party_cat_semistrong=="Independent"]


#ideology identification
fData <- fData[,Ideo_cat := ifelse(polviews%in%c("liberal","extremely liberal","slightly liberal"),"Liberal",
                                   ifelse(polviews%in%c("conservative","extremely conservative","slightly conservative"),"Conservative",
                                          ifelse(polviews=="moderate, middle of the road","Moderate",NA)))]

fData <- fData[,Lib := Ideo_cat=="Liberal"]
fData <- fData[,Cons := Ideo_cat=="Conservative"]
fData <- fData[,Mod := Ideo_cat=="Moderate"]


fData <- fData[,Ideo_cat_strong := ifelse(polviews%in%c("extremely liberal"),"Liberal",
                                          ifelse(polviews%in%c("extremely conservative"),"Conservative",
                                                 ifelse(polviews=="moderate, middle of the road","Moderate",NA)))]

fData <- fData[,Lib_strong := Ideo_cat_strong=="Liberal"]
fData <- fData[,Cons_strong := Ideo_cat_strong=="Conservative"]
fData <- fData[,Mod_strong := Ideo_cat_strong=="Moderate"]

fData <- fData[,Ideo_cat_semistrong := ifelse(polviews%in%c("liberal","extremely liberal"),"Liberal",
                                              ifelse(polviews%in%c("conservative","extremely conservative"),"Conservative",
                                                     ifelse(polviews=="moderate, middle of the road","Moderate",NA)))]

fData <- fData[,Lib_semistrong := Ideo_cat_semistrong=="Liberal"]
fData <- fData[,Cons_semistrong := Ideo_cat_semistrong=="Conservative"]
fData <- fData[,Mod_semistrong := Ideo_cat_semistrong=="Moderate"]

#get numerical party identity and ideology (for graphs)
fData <- fData[partyid=="strong republican",Party_id_numeric := 7]
fData <- fData[partyid=="not very strong republican",Party_id_numeric := 6]
fData <- fData[partyid=="independent, close to republican",Party_id_numeric := 5]
fData <- fData[partyid=="independent (neither, no response)",Party_id_numeric := 4]
fData <- fData[partyid=="independent, close to democrat",Party_id_numeric := 3]
fData <- fData[partyid=="not very strong democrat",Party_id_numeric := 2]
fData <- fData[partyid=="strong democrat",Party_id_numeric := 1]


#clean covariates ("" should be NA)
fData <- fData[partyid=="other party",partyid:=NA]
fData <- fData[partyid=="other party",Dem:=NA]
fData <- fData[partyid=="other party",reps:=NA]
# fData <- fData[agegroup=="",agegroup:=NA]
# fData <- fData[income5anes2=="",income5anes2:=NA]
fData <- fData[female=="",female:=NA]
fData <- fData[wrkstat=="",wrkstat:=NA]
fData <- fData[degree=="",degree:=NA]
fData <- fData[race=="",race:=NA]
fData <- fData[relig=="",relig:=NA]
fData <- fData[,Unemployed := wrkstat=="unemployed, laid off, looking for work"]
fData <- fData[,College := degree%in%c("bachelor's","gradaute")]
fData <- fData[,Black := race=="black"]
fData <- fData[degree=="less than high school",Degree_num := 1]
fData <- fData[degree=="high school",Degree_num := 2]
fData <- fData[degree=="associate/junior college",Degree_num := 3]
fData <- fData[degree=="bachelor's",Degree_num := 4]
fData <- fData[degree=="graduate",Degree_num := 5]
fData <- fData[,High_education := ifelse(Degree_num%in%4:5,"High education","Low education")]


#get party and name of president
fData <- fData[,Incumbent_party:=ifelse(year%in%c(1969:1976,1981:1992,2001:2008,2017:2020),"Republican","Democrat")]


fData <- fData[year%in%1969:1974,Incumbent_president:="Nixon"]
fData <- fData[year%in%1975:1976,Incumbent_president:="Ford"]
fData <- fData[year%in%1977:1980,Incumbent_president:="Carter"]
fData <- fData[year%in%1981:1988,Incumbent_president:="Reagan"]
fData <- fData[year%in%1989:1992,Incumbent_president:="Bush sr."]
fData <- fData[year%in%1993:2000,Incumbent_president:="Clinton"]
fData <- fData[year%in%2001:2008,Incumbent_president:="Bush jr."]
fData <- fData[year%in%2009:2016,Incumbent_president:="Obama"]
fData <- fData[year%in%2017:2020,Incumbent_president:="Trump"]
fData <- fData[year%in%2021:2024,Incumbent_president:="Biden"]

#create confidence variables
fData <- fData[,trustfed := ifelse(confed=="a great deal",1,
                                   ifelse(confed=="only some",0.5,
                                          ifelse(confed=="hardly any",0,NA)))]

fData <- fData[,trustlegisl := ifelse(conlegis=="a great deal",1,
                                      ifelse(conlegis=="only some",0.5,
                                             ifelse(conlegis=="hardly any",0,NA)))]



#Get 'own pres' variable
fData <- fData[,ownpres := ifelse((Rep==T&Incumbent_party=="Republican")|Dem==T&Incumbent_party=="Democrat",1,0)]
fData <- fData[,ownpresIdeo := ifelse((Cons==T&Incumbent_party=="Republican")|Lib==T&Incumbent_party=="Democrat",1,0)]

#get dummies for whether president is rep/dem
fData <- fData[,President_rep := as.numeric(Incumbent_party=="Republican")]
fData <- fData[,President_dem := as.numeric(Incumbent_party=="Democrat")]

fData <- fData[,Dem := as.numeric(Dem)]
fData <- fData[,Rep := as.numeric(Rep)]


#remove data with missing trust in federal gov
fData <- fData[!is.na(trustfed)]


#Majority party house and senate
fData <- fData[,Incumbent_house_party := ifelse(year%in%c(1995:2006,2011:2018),"Republican",
                                                ifelse(year%in%c(1973:1994,2007:2010,2019:2023),"Democrat",NA))]

fData <- fData[,Incumbent_senate_party := ifelse(year%in%c(1981:1986,1995:2000,2003:2006,2015:2020),"Republican",
                                                 ifelse(year%in%c(1973:1980,1987:1994,2001:2002,
                                                                  2007:2014,2021:2022),"Democrat",NA))]

#years of turnover elections (president)
vTurnover_elections <- c(1976,1980,1992,2000,2008,2016,2020) #years of turnover elections
fElections <- fData[,.(Incumbent_president = Incumbent_president[1],Incumbent_party = Incumbent_party[1]),by=year ] #all waves and incumbent presidents during that wave
vMinYear <- vMinYear_all[vMinYear_all%in%fElections$year]


#years of turnover elections house
vTurnover_elections_house <- c(1994,2006,2010,2018) #years of turnover elections
fElections_house <- fData[,.(Incumbent_house_party = Incumbent_house_party[1]),by=year ] #all waves and incumbent presidents during that wave
vMinYear_house <- vMinYear_all[vMinYear_all%in%fElections_house$year]

#years of turnover elections senate
vTurnover_elections_senate <- c(1980,1986,1994,2000,2002,2006,2014,2020) #years of turnover elections
fElections_senate <- fData[,.(Incumbent_senate_party = Incumbent_senate_party[1]),by=year ] #all waves and incumbent presidents during that wave
vMinYear_senate <- vMinYear_all[vMinYear_all%in%fElections_senate$year]



fData <- fData[,ownsenate := ifelse((Rep==T&Incumbent_senate_party=="Republican")|Dem==T&Incumbent_senate_party=="Democrat",1,0)]
fData <- fData[,ownhouse := ifelse((Rep==T&Incumbent_house_party=="Republican")|Dem==T&Incumbent_house_party=="Democrat",1,0)]


############ Figure 1: Confidence in federal government ##########
fGraph_main <- fData[!is.na(Party_cat),.(Av_conf = mean(trustfed,na.rm=T),
                                         Obs=.N),by=.(year,Party_cat)]
fGraph_main <- fGraph_main[!is.na(Av_conf)]


pdf("Graphs/Trust_federalgov.pdf")
ggplot(fGraph_main, aes(x=year,y=Av_conf,fill=Party_cat,col=Party_cat)) +
  geom_rect(data=fGraph_main[1,],aes(xmin = 1973.5, xmax = 1976.5, ymin = -Inf, ymax = Inf), fill = "red", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 1976.5, xmax = 1980.5, ymin = -Inf, ymax = Inf), fill = "blue", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 1980.5, xmax = 1992.5, ymin = -Inf, ymax = Inf), fill = "red", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 1992.5, xmax = 2000.5, ymin = -Inf, ymax = Inf), fill = "blue", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 2000.5, xmax = 2008.5, ymin = -Inf, ymax = Inf), fill = "red", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 2008.5, xmax = 2016.5, ymin = -Inf, ymax = Inf), fill = "blue", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 2016.5, xmax = 2020.5, ymin = -Inf, ymax = Inf), fill = "red", alpha = 0.2,color=NA)+
  geom_rect(data=fGraph_main[1,],aes(xmin = 2020.5, xmax = 2022, ymin = -Inf, ymax = Inf), fill = "blue", alpha = 0.2,color=NA)+
  geom_point() +
  geom_line() +
  # facet_wrap(~variable_name,scales="free_y") +
  theme_bw() +
  scale_color_manual(values=c("blue","green","red")) +
  scale_x_continuous(expand = c(0,0))+
  labs(color="",fill="") +
  xlab("Year")+
  scale_y_continuous(name="Confidence in federal government",labels=scales::percent_format(accuracy = 1))+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()






############ Figure 2, Panel A Evolution of president-in-power effect (turnover elections) ###############
lYearSens_allturnovers <- lapply(vTurnover_elections, function(x){ #run regression for each turnover election
  
  #whenever possible, select last two years before turnover eleciton and first two years after. Otherwise, one year after/before
  #Below is the selection for the relevant years
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  #get restricted data for turnover election
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  #run regression
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  #get all relevant info in table format
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$NewPres <- fElections[year>x]$Incumbent_president[1]
  fInfo_all$OldPres <- fElections[year<=x]$Incumbent_president[nrow(fElections[year<=x])]
  
  return(fInfo_all)
})

#combine all turnover election regression
fInfo_all_allturnovers <- do.call(rbind.data.frame, lYearSens_allturnovers)

#get elections in right order
fInfo_all_allturnovers <- as.data.table(fInfo_all_allturnovers)[,Turnover := paste0(OldPres, " - ",NewPres)]
colnames(fInfo_all_allturnovers) <- c("Estimate","SE","T_val","P_val","Effect","NewPres","OldPres","Election")
fInfo_all_allturnovers <- fInfo_all_allturnovers[,Election := factor(Election, levels=as.character(unique(fInfo_all_allturnovers$Election)))]
fInfo_all_allturnovers <- fInfo_all_allturnovers[,Election_num := as.numeric(Election)]

#prepare data frame for graph
fInfo_all_allturnovers <- as.data.table(fInfo_all_allturnovers)[,CI_lower := Estimate-1.96*SE]
fInfo_all_allturnovers <- as.data.table(fInfo_all_allturnovers)[,CI_upper := Estimate+1.96*SE]
fInfo_all_allturnovers <- fInfo_all_allturnovers[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

pdf("Graphs/Window_estimates_allturnovers.pdf")
ggplot(fInfo_all_allturnovers, aes(x=Election_num,y=Estimate,fill=Effect,col=Effect)) +
  geom_point() +
  # geom_smooth(se=F) +
  geom_line()+
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Institution) +
  scale_x_continuous(breaks=1:7,minor_breaks = 1:7, labels = as.character(unique(fInfo_all_allturnovers$Election)))+
  theme_bw() +
  xlab("Election")+
  labs(fill="",col="")+
  theme(text=element_text(size=18), axis.text.x=element_text(angle=90, vjust = 0.3),legend.position = "bottom")
dev.off()




############ Figure 2, Panel B Evolution of president-in-power effect (Shrinking window) ###############
lYearSens_decl <- lapply(vMinYear, function(x){ #run regression from x to the end of the sample
  fData_rest <- fData[year>=x] #get data from year x to the end of hte sample
  
  #run regression pip effect
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  #obtain all info for graphs
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  
  return(fInfo_all)
  
})

#combine all windows
fInfo_all_decl <- do.call(rbind.data.frame, lYearSens_decl)

#prepare data for graph
colnames(fInfo_all_decl) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")
fInfo_all_decl <- as.data.table(fInfo_all_decl)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decl <- as.data.table(fInfo_all_decl)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decl <- fInfo_all_decl[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]


pdf("Graphs/Window_estimates_decl.pdf")
ggplot(fInfo_all_decl, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()


############ Figure 2, Panel C Evolution of president-in-power effect (shifting window) ###############
lYearSens_shift <- lapply(vMinYear, function(x){ #run regression from x to x+iWindowlength
  fData_rest <- fData[year%in%x:(x+iWindowLength)] #get data from year x to x+iWindowlength
  
  #regression for pip effect on restricted data set
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  #get all relevant regression info
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  
  return(fInfo_all)
  
})

#combine all windows
fInfo_all_shift <- do.call(rbind.data.frame, lYearSens_shift)

#prepare data frame for graph
colnames(fInfo_all_shift) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")
fInfo_all_shift <- as.data.table(fInfo_all_shift)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shift <- as.data.table(fInfo_all_shift)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shift <- fInfo_all_shift[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

#make graph
pdf("Graphs/Window_estimates_shift.pdf")
ggplot(fInfo_all_shift, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()


########### In text: Test whether effect has changed differently over time for dems and reps. Last estimate compared to all others ##############
lYearSens_decl_tables <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  return(reg_1)
  
})

lWaldLastDeclWindow_asym <- lapply(1:length(lYearSens_decl_tables), function(x){
  reg1 <- lYearSens_decl_tables[[x]]
  reg2 <- lYearSens_decl_tables[[length(lYearSens_decl_tables)]]
  
  #main effect
  beta_main <- coef(summary(reg1))["ownpres", "Estimate"]
  var_beta_main <- coef(summary(reg1))["ownpres", "Std. Error"]^2
  
  beta_main_lastperiod <- coef(summary(reg2))["ownpres", "Estimate"]
  var_beta_main_lastperiod <- coef(summary(reg2))["ownpres", "Std. Error"]^2
  
  Wald_main<- ((beta_main - beta_main_lastperiod)^2) / (var_beta_main + var_beta_main_lastperiod)
  Pval_main <- 1 - pchisq(Wald_main, df = 1)
  
  
  
  
  
  #asymetry
  beta_asym <- coef(summary(reg1))["Rep:ownpres", "Estimate"]
  var_beta_asym <- coef(summary(reg1))["Rep:ownpres", "Std. Error"]^2
  
  beta_asym_lastperiod <- coef(summary(reg2))["Rep:ownpres", "Estimate"]
  var_beta_asym_lastperiod <- coef(summary(reg2))["Rep:ownpres", "Std. Error"]^2
  
  Wald_asym<- ((beta_asym - beta_asym_lastperiod)^2) / (var_beta_asym + var_beta_asym_lastperiod)
  Pval_asym <- 1 - pchisq(Wald_asym, df = 1)
  
  fInfo <- data.frame(Pval_main = Pval_main, Pval_asym = Pval_asym, Minyear = vMinYear[x], type ="Declining window")
  return(fInfo)
  
})

fWaldLastDeclWindow_asym <- do.call(rbind, lWaldLastDeclWindow_asym)



lYearSens_shift_tables <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  return(reg_1)
  
})

lWaldLastShiftWindow_asym <- lapply(1:length(lYearSens_shift_tables), function(x){
  reg1 <- lYearSens_shift_tables[[x]]
  reg2 <- lYearSens_shift_tables[[length(lYearSens_shift_tables)]]
  
  #main effect
  beta_main <- coef(summary(reg1))["ownpres", "Estimate"]
  var_beta_main <- coef(summary(reg1))["ownpres", "Std. Error"]^2
  
  beta_main_lastperiod <- coef(summary(reg2))["ownpres", "Estimate"]
  var_beta_main_lastperiod <- coef(summary(reg2))["ownpres", "Std. Error"]^2
  
  Wald_main<- ((beta_main - beta_main_lastperiod)^2) / (var_beta_main + var_beta_main_lastperiod)
  Pval_main <- 1 - pchisq(Wald_main, df = 1)
  
  #asymmetry
  beta_asym <- coef(summary(reg1))["Rep:ownpres", "Estimate"]
  var_beta_asym <- coef(summary(reg1))["Rep:ownpres", "Std. Error"]^2
  
  beta_asym_lastperiod <- coef(summary(reg2))["Rep:ownpres", "Estimate"]
  var_beta_asym_lastperiod <- coef(summary(reg2))["Rep:ownpres", "Std. Error"]^2
  
  Wald_asym<- ((beta_asym - beta_asym_lastperiod)^2) / (var_beta_asym + var_beta_asym_lastperiod)
  Pval_asym <- 1 - pchisq(Wald_asym, df = 1)
  
  fInfo <- data.frame(Pval_main = Pval_main, Pval_asym = Pval_asym, Minyear = vMinYear[x], type ="Shifting window")
  return(fInfo)
  
})

fWaldLastShiftWindow_asym <- do.call(rbind,lWaldLastShiftWindow_asym)


lYearSens_allturnovers_tables <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  return(reg_1)
  
})

lWaldLastTurnover_asym <- lapply(1:length(lYearSens_allturnovers_tables), function(x){
  reg1 <- lYearSens_allturnovers_tables[[x]]
  reg2 <- lYearSens_allturnovers_tables[[length(lYearSens_allturnovers_tables)]]
  
  #main effect
  beta_main <- coef(summary(reg1))["ownpres", "Estimate"]
  var_beta_main <- coef(summary(reg1))["ownpres", "Std. Error"]^2
  
  beta_main_lastperiod <- coef(summary(reg2))["ownpres", "Estimate"]
  var_beta_main_lastperiod <- coef(summary(reg2))["ownpres", "Std. Error"]^2
  
  Wald_main<- ((beta_main - beta_main_lastperiod)^2) / (var_beta_main + var_beta_main_lastperiod)
  Pval_main <- 1 - pchisq(Wald_main, df = 1)
  
  #asymmetry
  beta_asym <- coef(summary(reg1))["Rep:ownpres", "Estimate"]
  var_beta_asym <- coef(summary(reg1))["Rep:ownpres", "Std. Error"]^2
  
  beta_asym_lastperiod <- coef(summary(reg2))["Rep:ownpres", "Estimate"]
  var_beta_asym_lastperiod <- coef(summary(reg2))["Rep:ownpres", "Std. Error"]^2
  
  Wald_asym<- ((beta_asym - beta_asym_lastperiod)^2) / (var_beta_asym + var_beta_asym_lastperiod)
  Pval_asym <- 1 - pchisq(Wald_asym, df = 1)
  
  fInfo <- data.frame(Pval_main = Pval_main, Pval_asym = Pval_asym, Turnover = as.character(unique(fInfo_all_allturnovers$Election))[x], type ="Turnover elections")
  return(fInfo)
  
})

fWaldLastTurnover_asym <- do.call(rbind,lWaldLastTurnover_asym)



########### Figure 3: Strength of partisan identification ###############
fGraph_extremefraction_party <- fData[!is.na(partyid),.(Tot_dem = sum(partyid%in%c("independent, close to democrat", "not very strong democrat","strong democrat")),
                                                        Tot_rep = sum(partyid%in%c("independent, close to republican", "not very strong republican","strong republican")),
                                                        Tot_strong_dem = sum(partyid%in%c("strong democrat")),
                                                        Tot_strong_rep = sum(partyid%in%c("strong republican")),
                                                        Obs=.N),by=year]
fGraph_extremefraction_party <- fGraph_extremefraction_party[,Fraction_strong_dem := Tot_strong_dem/Tot_dem]
fGraph_extremefraction_party <- fGraph_extremefraction_party[,Fraction_strong_rep := Tot_strong_rep/Tot_rep]
fGraph_extremefraction_party <- fGraph_extremefraction_party[,c("year","Fraction_strong_dem","Fraction_strong_rep")]

fGraph_extremefraction_party_long <- reshape2::melt(fGraph_extremefraction_party,id.vars=c("year"))
fGraph_extremefraction_party_long <- as.data.table(fGraph_extremefraction_party_long)[,variable := ifelse(variable=="Fraction_strong_dem","Democrats","Republicans")]

pdf("Graphs/Fraction_extreme_party.pdf") 
ggplot(fGraph_extremefraction_party_long, aes(x=year,y=value,fill=variable,col=variable)) +
  geom_point() +
  geom_smooth(se=F) +
  theme_bw() +
  xlab("Year") +
  scale_y_continuous(name="Strong partisans (%)", labels=scales::percent_format(accuracy=1))+
  labs(fill="",col="")+
  scale_color_manual(values=c("blue","red")) +
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()





############# Table 1: Effect of alignment on trust in government by education leve ##################
reg_educ_combined_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[degree=="less than high school"]$wtssall, data=fData[degree=="less than high school"])
reg_educ_combined_2 <- felm(trustfed~Dem+Rep+ownpres+ownpres+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[degree=="high school"]$wtssall, data=fData[degree=="high school"])
reg_educ_combined_3 <- felm(trustfed~Dem+Rep+ownpres+ownpres+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[degree=="associate/junior college"]$wtssall, data=fData[degree=="associate/junior college"])
reg_educ_combined_4 <- felm(trustfed~Dem+Rep+ownpres+ownpres+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[degree=="bachelor's"]$wtssall, data=fData[degree=="bachelor's"])
reg_educ_combined_5 <- felm(trustfed~Dem+Rep+ownpres+ownpres+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[degree=="graduate"]$wtssall, data=fData[degree=="graduate"])


Regs_main_educ_combined <- stargazer(reg_educ_combined_1,reg_educ_combined_2,reg_educ_combined_3,reg_educ_combined_4,reg_educ_combined_5,
                                     summary=F,
                                     font.size = "scriptsize",
                                     title=paste0("Effect of alignment on presidential support by education level"),
                                     label="Results_main_educ_combined",
                                     table.placement = "H",
                                     covariate.labels = c("Democrat","Republican","Support pres."),
                                     # covariate.labels = c("Liberal","Conservative","Support pres.","Cons. x Support pres."),
                                     omit = c("realinc","age","female","educ","wrkstat","childs","realinc","year"),
                                     omit.stat = c("ser","f","rsq"),
                                     add.lines = list(c("Year fixed effects",rep("Yes",5)),
                                                      c("Socioeconomic controls",rep("Yes",5))),
                                     column.labels = c("less than high school","high school","associate/junior college","bachelor's","graduate"))




############# Figure 4: Difference in education level between Republicans and Democrats ############
fGraph_educ <- fData[!is.na(Party_cat),.(Av_college = mean(degree=="graduate",na.rm=T),
                                         Av_educ = mean(Degree_num,na.rm=T)),by=.(Party_cat,year)]
fGraph_educ_wide_college <- reshape2::dcast(fGraph_educ, year~Party_cat, value.var = "Av_college")
fGraph_educ_wide_college <- as.data.table(fGraph_educ_wide_college)[,Difference := Republican-Democrat]

fGraph_educ_wide_educ <- reshape2::dcast(fGraph_educ, year~Party_cat, value.var = "Av_educ")
fGraph_educ_wide_educ <- as.data.table(fGraph_educ_wide_educ)[,Difference := Republican-Democrat]


pdf("Graphs/Education_level_difference.pdf")
ggplot(fGraph_educ_wide_educ, aes(x=year, y=Difference)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  geom_hline(yintercept = 0)+
  labs(fill="",col="") +
  ylab("Difference in education level Reps and Dems") +
  xlab("Year")+
  theme(text=element_text(size=15), legend.position = "bottom")
dev.off()


########### Table A1: Summary statistics ##############
fData_sumstats <- fData[!is.na(wtssall)&!is.na(trustfed)&Party_cat%in%c("Democrat","Republican","Independent")]
vCovariates <- c("trustfed", "age","Unemployed","realinc","College","Black")


fSumStats_party <- fData_sumstats[,lapply(.SD, function(x) sprintf("%.2f",mean(x,na.rm=T))),by=.(Party_cat,Incumbent_party),.SDcols=c(vCovariates)]
fSumStats_party <- fSumStats_party[order(Party_cat)]

fSumStats_party_t <- as.data.frame(t(fSumStats_party))
colnames(fSumStats_party_t) <- paste0(as.matrix(fSumStats_party_t[1,1:6]))
fSumStats_party_t <- fSumStats_party_t[-1,]
rownames(fSumStats_party_t) <- c("Party president","Confidence Federal Government","Age","Unemployed","Income","College","Black")
colnames(fSumStats_party_t) <- c("Dem_rep","Dem_dem","Ind_rep","Ind_dem", "Rep_rep","Rep_dem")

fSumStats_party_t_tex <- stargazer(fSumStats_party_t, 
                                   summary = F, 
                                   title="Summary statistics",
                                   font.size = "scriptsize",
                                   digits=2,
                                   table.placement = "H")
write.table(fSumStats_party_t_tex,file=paste0("Tables/Summary_stats_party)asymmetry.tex"),row.names= FALSE,na="", quote = FALSE, col.names = FALSE)



########### Figure A1 and Table A2: ############
#Create Figure A1
lYearSens_allturnovers_pretrends <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  #make year a factor, and set the year before the election as base
  fData_rest <- fData_rest[,Year_factor := as.factor(as.numeric(year))]
  fData_rest <- fData_rest[,Year_factor := relevel(Year_factor,ref=ifelse(length(vYears_before)==2,2,1))]
  
  fData_rest <- fData_rest[,Dem := as.numeric(Dem)]
  fData_rest <- fData_rest[,Rep := as.numeric(Rep)]
  
  # reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  #estimate effect per year
  reg_2 <- felm(trustfed~Dem+Rep+Dem:Year_factor+Rep:Year_factor+realinc|Year_factor+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  fCoefs <- as.data.frame(coef(summary(reg_2)))
  fCoefs <- fCoefs[grepl("Year_factor",rownames(fCoefs)),]
  fCoefs$Year <- gsub("Dem:Year_factor|Rep:Year_factor","",rownames(fCoefs))
  fCoefs$Party <- substr(rownames(fCoefs),1,3)
  fCoefs <- as.data.table(fCoefs)
  fCoefs <- fCoefs[,Before := Year%in%vYears_before]
  fCoefs <- fCoefs[Before==T,Year_num := -2, by=Party]
  fCoefs <- fCoefs[Before==F,Year_num := 0:(.N-1), by=Party]
  fCoefs <- fCoefs[,CI_lower := Estimate-1.96*`Std. Error`]
  fCoefs <- fCoefs[,CI_upper := Estimate+1.96*`Std. Error`]
  
  #create 0s for -1 (reference)
  fCoefs_null <- data.frame(Estimate = c(0,0),
                            `Std. Error` = c(0,0),
                            `t value` = NA,
                            `Pr(>|t)`=NA,
                            Year=ifelse(length(vYears_before)==2,vYears_before[2],vYears_before[1]),
                            Party=c("Dem","Rep"),
                            Before=NA,
                            Year_num=-1,
                            CI_lower=0,
                            CI_upper=0)
  fCoefs_full <- rbind(fCoefs,fCoefs_null,use.names=F)
  fCoefs_full <- fCoefs_full[order(Party,Year_num)]
  
  fCoefs_full$NewPres <- fElections[year>x]$Incumbent_president[1]
  fCoefs_full$OldPres <- fElections[year<=x]$Incumbent_president[nrow(fElections[year<=x])]
  fCoefs_full <- fCoefs_full[,Turnover := paste0(OldPres, " - ",NewPres)]
  
  return(fCoefs_full)
  
})

fYearSens_allturnovers_pretrends <- do.call(rbind,lYearSens_allturnovers_pretrends)

fYearSens_allturnovers_pretrends <- fYearSens_allturnovers_pretrends[,Turnover := factor(Turnover, levels = unique(Turnover))]

#remove Trump - biden because there is only one pre-period (which is the reference)
fYearSens_allturnovers_pretrends <- fYearSens_allturnovers_pretrends[Turnover!="Trump - Biden"]

pdf("Graphs/Graph_did_pretrends_allturnovers.pdf")
ggplot(fYearSens_allturnovers_pretrends, aes(x=Year_num,y=Estimate,fill=Party,col=Party)) +
  geom_point() +
  geom_errorbar(aes(ymin=CI_lower,ymax=CI_upper)) +
  facet_wrap(~Turnover)+
  theme_bw() +
  geom_hline(yintercept = 0) +
  xlab("Years to election")
dev.off()



#create Table A2
lYearSens_allturnovers_pretrends_regs <- lapply(vTurnover_elections[vTurnover_elections!=2020], function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  print(vYears_before)
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  #make year a factor, and set the year before the election as base
  
  fData_rest <- fData_rest[order(year)]
  fData_rest <- fData_rest[,Year_factor := as.numeric(as.factor(year))-3]
  fData_rest <- fData_rest[,Year_factor := relevel(as.factor(Year_factor),ref=ifelse(length(vYears_before)==2,2,1))]
  
  fData_rest <- fData_rest[,Dem := as.numeric(Dem)]
  fData_rest <- fData_rest[,Rep := as.numeric(Rep)]
  
  # reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  #estimate effect per year
  reg_2 <- felm(trustfed~Dem+Rep+Dem:Year_factor+Rep:Year_factor+realinc|Year_factor+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  
  return(reg_2)
  
})


stargazer(lYearSens_allturnovers_pretrends_regs,
          summary=F,
          font.size = "tiny",
          title="President-in-power effect for per turnover election",
          label="Results_allturnovers_table_regs",
          column.labels = as.character(unique(fYearSens_allturnovers_pretrends$Turnover)),
          model.numbers = F,
          dep.var.caption = "",
          dep.var.labels="",
          # covariate.labels = c("Education"),
          covariate.labels = c("Democrat","Republican",
                               "Democrat x Year = -2",
                               "Democrat x Year = 0",
                               "Democrat x Year = 1",
                               "Republican x Year = -2",
                               "Republican x Year = 0",
                               "Republican x Year = 1"),
          omit = c("realinc","age","Female","educ","wrkstatoyed","childs","Married_bin","year"),
          omit.stat = c("ser","f","rsq"),
          add.lines = list(c("Year fixed effects",rep("Yes",length(lYearSens_allturnovers_pretrends_regs))),
                           c("Socioeconomic controls",rep("Yes",length(lYearSens_allturnovers_pretrends_regs)))))





########### Figure A2: Evolution of president-in-power effect, event study ##########

#Panel A
lYearSens_allturnoverseventstudy <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc+year+year:Rep|age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$NewPres <- fElections[year>x]$Incumbent_president[1]
  fInfo_all$OldPres <- fElections[year<=x]$Incumbent_president[nrow(fElections[year<=x])]
  return(fInfo_all)
  
})

fInfo_all_allturnoverseventstudy <- do.call(rbind.data.frame, lYearSens_allturnoverseventstudy)
fInfo_all_allturnoverseventstudy <- as.data.table(fInfo_all_allturnoverseventstudy)[,Turnover := paste0(OldPres, " - ",NewPres)]

colnames(fInfo_all_allturnoverseventstudy) <- c("Estimate","SE","T_val","P_val","Effect","NewPres","OldPres","Election")


#get right order
fInfo_all_allturnoverseventstudy <- fInfo_all_allturnoverseventstudy[,Election := factor(Election, levels=as.character(unique(fInfo_all_allturnoverseventstudy$Election)))]
fInfo_all_allturnoverseventstudy <- fInfo_all_allturnoverseventstudy[,Election_num := as.numeric(Election)]

fInfo_all_allturnoverseventstudy <- as.data.table(fInfo_all_allturnoverseventstudy)[,CI_lower := Estimate-1.96*SE]
fInfo_all_allturnoverseventstudy <- as.data.table(fInfo_all_allturnoverseventstudy)[,CI_upper := Estimate+1.96*SE]
fInfo_all_allturnoverseventstudy <- fInfo_all_allturnoverseventstudy[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

pdf("Graphs/Window_estimates_allturnoverseventstudy.pdf")
ggplot(fInfo_all_allturnoverseventstudy, aes(x=Election_num,y=Estimate,fill=Effect,col=Effect)) +
  geom_point() +
  # geom_smooth(se=F) +
  geom_line()+
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Institution) +
  scale_x_continuous(breaks=1:7,minor_breaks = 1:7, labels = as.character(unique(fInfo_all_allturnoverseventstudy$Election)))+
  theme_bw() +
  xlab("Election")+
  labs(fill="",col="")+
  theme(text=element_text(size=18), axis.text.x=element_text(angle=90, vjust = 0.3),legend.position = "bottom")
dev.off()



#Panel B
lYearSens_decleventstudy <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  fData_rest <- fData_rest[Party_cat%in%c("Republican","Democrat")]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc+year+year:Rep|age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_decleventstudy <- do.call(rbind.data.frame, lYearSens_decleventstudy)

colnames(fInfo_all_decleventstudy) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")

fInfo_all_decleventstudy <- as.data.table(fInfo_all_decleventstudy)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decleventstudy <- as.data.table(fInfo_all_decleventstudy)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decleventstudy <- fInfo_all_decleventstudy[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]


pdf("Graphs/Window_estimates_decleventstudy.pdf")
ggplot(fInfo_all_decleventstudy, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()


#Panel C
lYearSens_shifteventstudy <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc+year+year:Rep|age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_shifteventstudy <- do.call(rbind.data.frame, lYearSens_shifteventstudy)

colnames(fInfo_all_shifteventstudy) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")

fInfo_all_shifteventstudy <- as.data.table(fInfo_all_shifteventstudy)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shifteventstudy <- as.data.table(fInfo_all_shifteventstudy)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shifteventstudy <- fInfo_all_shifteventstudy[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]


pdf("Graphs/Window_estimates_shifteventstudy.pdf")
ggplot(fInfo_all_shifteventstudy, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()



######### Table A3: President-in-power effect for conservative Democrats and liberal Republicans #############
fParty_vs_ideo <- fData[!is.na(Cons)&!is.na(Lib)&!is.na(Rep)&!is.na(Dem),.(Av_dem = mean(Dem_semistrong,na.rm=T),
                                                                           Av_rep = mean(Rep_semistrong,na.rm=T),
                                                                           Obs=.N),
                        by=.(Ideo_cat)]


fData_consdem_librep <- fData[(Cons==T&Dem==T)|(Lib==T&Rep==T)|Ind==T] #data for conservative Democrats and liberal Republicans

reg_consdem_librep_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_consdem_librep$wtssall, data=fData_consdem_librep)


Regs_consdem_librep <- stargazer(reg_consdem_librep_1,
                                 summary=F,
                                 font.size = "tiny",
                                 title=paste0("President-in-power effect for conservative Democrats and liberal Republicans"),
                                 label="Results_consdem_librep",
                                 table.placement = "H",
                                 # covariate.labels = c("Education"),
                                 covariate.labels = c("Democrat","Republican","President-in-power effect","Rep. x President-in-power effect"),
                                 omit = c("realinc","age","Female","educ","wrkstatoyed","childs","Married_bin","year"),
                                 omit.stat = c("ser","f","rsq"),
                                 add.lines = list(c("Year fixed effects","Yes"),
                                                  c("Socioeconomic controls","Yes")))




########### Figure A3: Perceived ideological distance to parties ###########
#import data anes
fData_anes <- fread("ANES_clean.csv")

fGraph_perception_distance <- fData_anes[!is.na(Party_cat)&Party_cat!="",.(Perceived_distance_dem = mean(Perceived_distance_dem,na.rm=T),
                                                                      Perceived_distance_rep = mean(Perceived_distance_rep,na.rm=T)),by=.(Year,Party_cat)]
fGraph_perception_distance <- fGraph_perception_distance[!is.na(Perceived_distance_dem)]
fGraph_perception_distance_long <- reshape2::melt(fGraph_perception_distance, id.vars=c("Year","Party_cat"))
fGraph_perception_distance_long <- as.data.table(fGraph_perception_distance_long)[,Perception_party := ifelse(variable=="Perceived_distance_dem","Perceived distance to Dems","Perceived distance to Reps")]

pdf("Graphs/ANES_perception_ideology_distance.pdf")
ggplot(fGraph_perception_distance_long,aes(x=Year, y=value,fill=Party_cat,col=Party_cat)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  facet_wrap(~Perception_party) +
  scale_color_manual(values=c("blue","lightgreen","red")) +
  theme(text=element_text(size=15),legend.position = "bottom") +
  ylab("Distance") +
  labs(fill="",col="")
dev.off()


########### Table A4: President-in-power effect per turnover election ############
lYearSens_allturnovers_tables <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  return(reg_1)
  
})

stargazer(lYearSens_allturnovers_tables,
          summary=F,
          font.size = "tiny",
          title="President-in-power effect for per turnover election",
          label="Results_allturnovers_table",
          column.labels = as.character(unique(fInfo_all_allturnovers$Turnover)),
          model.numbers = F,
          dep.var.caption = "",
          dep.var.labels="",
          # covariate.labels = c("Education"),
          covariate.labels = c("Democrat","Republican","President-in-power effect","Rep. x President-in-power effect"),
          omit = c("realinc","age","Female","educ","wrkstatoyed","childs","Married_bin","year"),
          omit.stat = c("ser","f","rsq"),
          add.lines = list(c("Year fixed effects",rep("Yes",length(lYearSens_allturnovers_tables))),
                           c("Socioeconomic controls",rep("Yes",length(lYearSens_allturnovers_tables)))))




########## Table A5: President-in-power effect by the strength of partisan identification ###########

#estimate effect for independent, close to party, not very strong party, and strong party
reg_partyid_weak <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[!partyid%in%c("strong democrat","strong republican","not very strong democrat","not very strong republican")]$wtssall,data=fData[!partyid%in%c("strong democrat","strong republican","not very strong democrat","not very strong republican")])
reg_partyid_medium <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[!partyid%in%c("strong democrat","strong republican","independent, close to democrat","independent, close to republican")]$wtssall, data=fData[!partyid%in%c("strong democrat","strong republican","independent, close to democrat","independent, close to republican")])
reg_partyid_strong <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData[!partyid%in%c("independent, close to democrat","independent, close to republican","not very strong democrat","not very strong republican")]$wtssall, data= fData[!partyid%in%c("independent, close to democrat","independent, close to republican","not very strong democrat","not very strong republican")])


Regs_partyid <- stargazer(reg_partyid_weak, reg_partyid_medium, reg_partyid_strong,
                          summary=F,
                          font.size = "scriptsize",
                          title=paste0("President-in-power effect by strength of partisan identification"),
                          label="Results_partyid_table",
                          table.placement = "H",
                          model.numbers = F,
                          dep.var.caption = "",
                          column.labels = c("Weak identifiers","Medium identifiers","Strong identifiers"),
                          covariate.labels = c("Democrat","Republican","President-in-power effect","Rep. x President-in-power effect"),
                          omit = c("realinc","age","Female","educ","wrkstatoyed","childs","Married_bin","year"),
                          omit.stat = c("ser","f","rsq"),
                          add.lines = list(c("Year fixed effects","Yes","Yes","Yes"),
                                           c("Socioeconomic controls","Yes","Yes","Yes")))




############ Figure A4: Extreme partisans ###############
#Panel A
#per-turnover window routine for extreme partisans
lYearSens_allturnovers_extremes <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  reg_1 <- felm(trustfed~Dem_strong+Rep_strong+ownpres+ownpres:Rep_strong+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep_strongTRUE:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$NewPres <- fElections[year>x]$Incumbent_president[1]
  fInfo_all$OldPres <- fElections[year<=x]$Incumbent_president[nrow(fElections[year<=x])]
  return(fInfo_all)
  
})

fInfo_all_allturnovers_extremes <- do.call(rbind.data.frame, lYearSens_allturnovers_extremes)
fInfo_all_allturnovers_extremes <- as.data.table(fInfo_all_allturnovers_extremes)[,Turnover := paste0(OldPres, " - ",NewPres)]

colnames(fInfo_all_allturnovers_extremes) <- c("Estimate","SE","T_val","P_val","Effect","NewPres","OldPres","Election")


#get right order
fInfo_all_allturnovers_extremes <- fInfo_all_allturnovers_extremes[,Election := factor(Election, levels=as.character(unique(fInfo_all_allturnovers_extremes$Election)))]
fInfo_all_allturnovers_extremes <- fInfo_all_allturnovers_extremes[,Election_num := as.numeric(Election)]

fInfo_all_allturnovers_extremes <- as.data.table(fInfo_all_allturnovers_extremes)[,CI_lower := Estimate-1.96*SE]
fInfo_all_allturnovers_extremes <- as.data.table(fInfo_all_allturnovers_extremes)[,CI_upper := Estimate+1.96*SE]
fInfo_all_allturnovers_extremes <- fInfo_all_allturnovers_extremes[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

pdf("Graphs/Window_estimates_allturnovers_extremes.pdf")
ggplot(fInfo_all_allturnovers_extremes, aes(x=Election_num,y=Estimate,fill=Effect,col=Effect)) +
  geom_point() +
  # geom_smooth(se=F) +
  geom_line()+
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Institution) +
  scale_x_continuous(breaks=1:7,minor_breaks = 1:7, labels = as.character(unique(fInfo_all_allturnovers_extremes$Election)))+
  theme_bw() +
  xlab("Election")+
  labs(fill="",col="")+
  theme(text=element_text(size=18), axis.text.x=element_text(angle=90, vjust = 0.3),legend.position = "bottom")
dev.off()


#Panel B
#declining window routine for extreme partisans
lYearSens_decl_extremes <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  
  reg_1 <- felm(trustfed~Dem_strong+Rep_strong+ownpres+ownpres:Rep_strong+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep_strongTRUE:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_decl_extremes <- do.call(rbind.data.frame, lYearSens_decl_extremes)

colnames(fInfo_all_decl_extremes) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")

fInfo_all_decl_extremes <- as.data.table(fInfo_all_decl_extremes)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decl_extremes <- as.data.table(fInfo_all_decl_extremes)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decl_extremes <- fInfo_all_decl_extremes[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]


pdf("Graphs/Window_estimates_decl_extremes.pdf")
ggplot(fInfo_all_decl_extremes, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()

#Panel C
#Shfiting window routine for extreme partisans
lYearSens_shift_extremes <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  reg_1 <- felm(trustfed~Dem_strong+Rep_strong+ownpres+ownpres:Rep_strong+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep_strongTRUE:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_shift_extremes <- do.call(rbind.data.frame, lYearSens_shift_extremes)

colnames(fInfo_all_shift_extremes) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")

fInfo_all_shift_extremes <- as.data.table(fInfo_all_shift_extremes)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shift_extremes <- as.data.table(fInfo_all_shift_extremes)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shift_extremes <- fInfo_all_shift_extremes[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]


pdf("Graphs/Window_estimates_shift_extremes.pdf")
ggplot(fInfo_all_shift_extremes, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()










######## Figure A5: Alternative partisanship classification ##########
#Panel A
#Per-turnover election routine for alternative partisan classification
lYearSens_allturnovers_altclass <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  reg_1 <- felm(trustfed~Dem_semistrong+Rep_semistrong+ownpres+ownpres:Rep_semistrong+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep_semistrongTRUE:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$NewPres <- fElections[year>x]$Incumbent_president[1]
  fInfo_all$OldPres <- fElections[year<=x]$Incumbent_president[nrow(fElections[year<=x])]
  return(fInfo_all)
  
})

fInfo_all_allturnovers_altclass <- do.call(rbind.data.frame, lYearSens_allturnovers_altclass)
fInfo_all_allturnovers_altclass <- as.data.table(fInfo_all_allturnovers_altclass)[,Turnover := paste0(OldPres, " - ",NewPres)]

colnames(fInfo_all_allturnovers_altclass) <- c("Estimate","SE","T_val","P_val","Effect","NewPres","OldPres","Election")


#get right order
fInfo_all_allturnovers_altclass <- fInfo_all_allturnovers_altclass[,Election := factor(Election, levels=as.character(unique(fInfo_all_allturnovers_altclass$Election)))]
fInfo_all_allturnovers_altclass <- fInfo_all_allturnovers_altclass[,Election_num := as.numeric(Election)]

fInfo_all_allturnovers_altclass <- as.data.table(fInfo_all_allturnovers_altclass)[,CI_lower := Estimate-1.96*SE]
fInfo_all_allturnovers_altclass <- as.data.table(fInfo_all_allturnovers_altclass)[,CI_upper := Estimate+1.96*SE]
fInfo_all_allturnovers_altclass <- fInfo_all_allturnovers_altclass[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

pdf("Graphs/Window_estimates_allturnovers_altclass.pdf")
ggplot(fInfo_all_allturnovers_altclass, aes(x=Election_num,y=Estimate,fill=Effect,col=Effect)) +
  geom_point() +
  # geom_smooth(se=F) +
  geom_line()+
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Institution) +
  scale_x_continuous(breaks=1:7,minor_breaks = 1:7, labels = as.character(unique(fInfo_all_allturnovers_altclass$Election)))+
  theme_bw() +
  xlab("Election")+
  labs(fill="",col="")+
  theme(text=element_text(size=18), axis.text.x=element_text(angle=90, vjust = 0.3),legend.position = "bottom")
dev.off()


#Panel B
#Declining window routine for alternative partisan classification
lYearSens_decl_altclass <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  
  reg_1 <- felm(trustfed~Dem_semistrong+Rep_semistrong+ownpres+ownpres:Rep_semistrong+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep_semistrongTRUE:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_decl_altclass <- do.call(rbind.data.frame, lYearSens_decl_altclass)

colnames(fInfo_all_decl_altclass) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")

fInfo_all_decl_altclass <- as.data.table(fInfo_all_decl_altclass)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decl_altclass <- as.data.table(fInfo_all_decl_altclass)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decl_altclass <- fInfo_all_decl_altclass[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

pdf("Graphs/Window_estimates_decl_altclass.pdf")
ggplot(fInfo_all_decl_altclass, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()

#Panel C
#Shifting window routine for alternative partisan classification
lYearSens_shift_altclass <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  reg_1 <- felm(trustfed~Dem_semistrong+Rep_semistrong+ownpres+ownpres:Rep_semistrong+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep_semistrongTRUE:ownpres",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
  fInfo_all$Effect <- c("Main effect","Asymmetry")
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_shift_altclass <- do.call(rbind.data.frame, lYearSens_shift_altclass)

colnames(fInfo_all_shift_altclass) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")

fInfo_all_shift_altclass <- as.data.table(fInfo_all_shift_altclass)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shift_altclass <- as.data.table(fInfo_all_shift_altclass)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shift_altclass <- fInfo_all_shift_altclass[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]

pdf("Graphs/Window_estimates_shift_altclass.pdf")
ggplot(fInfo_all_shift_altclass, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  # facet_wrap(~Effect) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()










########### Figure A6: Evolution of president-in-power effect per education level, turnover elections###############
vDegrees <- unique(fData[!is.na(degree)]$degree)


#Per-turnover election routine, for each education group seperately
lYearSensallturnovers_educ_cat <- lapply(vTurnover_elections, function(x){
  if(x%in%c(2016)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }else if(x%in%c(2020)){
    vYears_after <- fElections[year>x]$year[1]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])):nrow(fElections[year<=x])]
  }else{
    vYears_after <- fElections[year>x]$year[1:2]
    vYears_before <- fElections[year<=x]$year[(nrow(fElections[year<=x])-1):nrow(fElections[year<=x])]
  }
  fData_rest <- fData[year%in%c(vYears_after,vYears_before)]
  
  lInfo_all <- lapply(vDegrees, function(sVar){
    fData_rest <- fData_rest[degree==sVar]
    reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
    fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
    fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
    
    fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
    fInfo_all$Effect <- c("Main effect","Asymmetry")
    fInfo_all$NewPres <- fElections[year>x]$Incumbent_president[1]
    fInfo_all$OldPres <- fElections[year<=x]$Incumbent_president[nrow(fElections[year<=x])]
    fInfo_all$Degree<- sVar
    return(fInfo_all)
  })
  
  fInfo_all_combined <- do.call(rbind,lInfo_all)
  
  return(fInfo_all_combined)
  
})

fInfo_all_allturnovers_educ_cat <- do.call(rbind.data.frame, lYearSensallturnovers_educ_cat)
fInfo_all_allturnovers_educ_cat <- as.data.table(fInfo_all_allturnovers_educ_cat)[,Turnover := paste0(OldPres, " - ",NewPres)]

colnames(fInfo_all_allturnovers_educ_cat) <- c("Estimate","SE","T_val","P_val","Effect","NewPres","OldPres","Degree","Election")

fInfo_all_allturnovers_educ_cat <- fInfo_all_allturnovers_educ_cat[,Election := factor(Election, levels=as.character(unique(fInfo_all_allturnovers_educ_cat$Election)))]
fInfo_all_allturnovers_educ_cat <- fInfo_all_allturnovers_educ_cat[,Election_num := as.numeric(Election)]

fInfo_all_allturnovers_educ_cat <- as.data.table(fInfo_all_allturnovers_educ_cat)[,CI_lower := Estimate-1.96*SE]
fInfo_all_allturnovers_educ_cat <- as.data.table(fInfo_all_allturnovers_educ_cat)[,CI_upper := Estimate+1.96*SE]
fInfo_all_allturnovers_educ_cat <- fInfo_all_allturnovers_educ_cat[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]
fInfo_all_allturnovers_educ_cat <- fInfo_all_allturnovers_educ_cat[,Degree := factor(Degree, levels=c("less than high school","high school","associate/junior college","bachelor's","graduate"))]


pdf("Graphs/Window_estimatesallturnovers_educ_cat.pdf")
ggplot(fInfo_all_allturnovers_educ_cat, aes(x=Election_num,y=Estimate,fill=Effect,col=Effect)) +
  geom_point() +
  # geom_smooth(se=F) +
  geom_line()+
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Degree) +
  scale_x_continuous(breaks=1:7,minor_breaks = 1:7, labels = as.character(unique(fInfo_all_allturnovers_educ_cat$Election)))+
  theme_bw() +
  xlab("Election")+
  labs(fill="",col="")+
  theme(text=element_text(size=18), axis.text.x=element_text(angle=90, vjust = 0.3),legend.position = "bottom")
dev.off()



############## Figure A7: Evolution of president-in-power effect per education level, shrinking window ###############
vDegrees <- unique(fData[!is.na(degree)]$degree)

#Declining window routine, for each education group seperately
lYearSens_decl_educ_cat <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  
  lInfo_all <- lapply(vDegrees, function(sVar){
    fData_rest <- fData_rest[degree==sVar]
    reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
    fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
    fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
    
    fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
    fInfo_all$Effect <- c("Main effect","Asymmetry")
    fInfo_all$MinYear <- x
    fInfo_all$Degree<- sVar
    return(fInfo_all)
  })
  
  fInfo_all_combined <- do.call(rbind,lInfo_all)
  
  return(fInfo_all_combined)
  
})

fInfo_all_decl_educ_cat <- do.call(rbind.data.frame, lYearSens_decl_educ_cat)

colnames(fInfo_all_decl_educ_cat) <- c("Estimate","SE","T_val","P_val","Effect","MinYear","Degree")

fInfo_all_decl_educ_cat <- as.data.table(fInfo_all_decl_educ_cat)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decl_educ_cat <- as.data.table(fInfo_all_decl_educ_cat)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decl_educ_cat <- fInfo_all_decl_educ_cat[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]
fInfo_all_decl_educ_cat <- fInfo_all_decl_educ_cat[,Degree := factor(Degree, levels=c("less than high school","high school","associate/junior college","bachelor's","graduate"))]

pdf("Graphs/Window_estimates_decl_educ_cat.pdf")
ggplot(fInfo_all_decl_educ_cat, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Degree) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()





########### Figure A8: Evolution of president-in-power effect per education level, shifting window ###########
vDegrees <- unique(fData[!is.na(degree)]$degree)


#Shifting window routine, for each education group seperately
lYearSens_shift_educ_cat <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  lInfo_all <- lapply(vDegrees, function(sVar){
    fData_rest <- fData_rest[degree==sVar]
    reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
    fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
    fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
    
    fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
    fInfo_all$Effect <- c("Main effect","Asymmetry")
    fInfo_all$MinYear <- x
    fInfo_all$Degree<- sVar
    return(fInfo_all)
  })
  
  fInfo_all_combined <- do.call(rbind,lInfo_all)
  
  return(fInfo_all_combined)
  
})

fInfo_all_shift_educ_cat <- do.call(rbind.data.frame, lYearSens_shift_educ_cat)

colnames(fInfo_all_shift_educ_cat) <- c("Estimate","SE","T_val","P_val","Effect","MinYear","Degree")

fInfo_all_shift_educ_cat <- as.data.table(fInfo_all_shift_educ_cat)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shift_educ_cat <- as.data.table(fInfo_all_shift_educ_cat)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shift_educ_cat <- fInfo_all_shift_educ_cat[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]
fInfo_all_shift_educ_cat <- fInfo_all_shift_educ_cat[,Degree := factor(Degree, levels=c("less than high school","high school","associate/junior college","bachelor's","graduate"))]

pdf("Graphs/Window_estimates_shift_educ_cat.pdf")
ggplot(fInfo_all_shift_educ_cat, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Degree) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()




############# Figure A9: Evolution of president/congress-in-power effect #############

#Panel A
#alignment on conf in fed gov, declining window
lYearSens_decl_congress_conffed <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+ownhouse+ownhouse:Rep+ownsenate+ownsenate:Rep++realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est_pres <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_levelpip_est_senate <- as.data.frame(summary(reg_1)$coefficients)["ownsenate",]
  fInfo_1_levelpip_est_house <- as.data.frame(summary(reg_1)$coefficients)["ownhouse",]
  fInfo_1_asymm_est_pres <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  fInfo_1_asymm_est_senate <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownsenate",]
  fInfo_1_asymm_est_house <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownhouse",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est_pres,fInfo_1_levelpip_est_senate,fInfo_1_levelpip_est_house,fInfo_1_asymm_est_pres,fInfo_1_asymm_est_senate,fInfo_1_asymm_est_house))
  fInfo_all$Effect <- c("Main effect President","Main effect Senate","Main effect House","Asymmetry President","Asymmetry Senate","Asymmetry House")
  fInfo_all$Government <- ifelse(grepl("Senate",fInfo_all$Effect),"Senate",
                                 ifelse(grepl("House",fInfo_all$Effect),"House","President"))
  fInfo_all$Effect <- gsub(" Senate| House| President","",fInfo_all$Effect)
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_decl_congress_conffed <- do.call(rbind.data.frame, lYearSens_decl_congress_conffed)

colnames(fInfo_all_decl_congress_conffed) <- c("Estimate","SE","T_val","P_val","Effect","Government", "MinYear")

fInfo_all_decl_congress_conffed <- as.data.table(fInfo_all_decl_congress_conffed)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decl_congress_conffed <- as.data.table(fInfo_all_decl_congress_conffed)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decl_congress_conffed <- fInfo_all_decl_congress_conffed[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]



pdf("Graphs/Window_estimates_decl_congress_conffed.pdf")
ggplot(fInfo_all_decl_congress_conffed, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Government) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()


#Panel B
#alignment on conf in fed gov, shifting window
lYearSens_shift_congress_conffed <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+ownhouse+ownhouse:Rep+ownsenate+ownsenate:Rep++realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est_pres <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_levelpip_est_senate <- as.data.frame(summary(reg_1)$coefficients)["ownsenate",]
  fInfo_1_levelpip_est_house <- as.data.frame(summary(reg_1)$coefficients)["ownhouse",]
  fInfo_1_asymm_est_pres <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  fInfo_1_asymm_est_senate <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownsenate",]
  fInfo_1_asymm_est_house <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownhouse",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est_pres,fInfo_1_levelpip_est_senate,fInfo_1_levelpip_est_house,fInfo_1_asymm_est_pres,fInfo_1_asymm_est_senate,fInfo_1_asymm_est_house))
  fInfo_all$Effect <- c("Main effect President","Main effect Senate","Main effect House","Asymmetry President","Asymmetry Senate","Asymmetry House")
  fInfo_all$Government <- ifelse(grepl("Senate",fInfo_all$Effect),"Senate",
                                 ifelse(grepl("House",fInfo_all$Effect),"House","President"))
  fInfo_all$Effect <- gsub(" Senate| House| President","",fInfo_all$Effect)
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_shift_congress_conffed <- do.call(rbind.data.frame, lYearSens_shift_congress_conffed)

colnames(fInfo_all_shift_congress_conffed) <- c("Estimate","SE","T_val","P_val","Effect","Government", "MinYear")

fInfo_all_shift_congress_conffed <- as.data.table(fInfo_all_shift_congress_conffed)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shift_congress_conffed <- as.data.table(fInfo_all_shift_congress_conffed)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shift_congress_conffed <- fInfo_all_shift_congress_conffed[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]



pdf("Graphs/Window_estimates_shift_congress_conffed_conffed.pdf")
ggplot(fInfo_all_shift_congress_conffed, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Government) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()


#Panel C
#alignment on conf in congress, declining window
lYearSens_decl_congress_conflegisl <- lapply(vMinYear, function(x){
  fData_rest <- fData[year>=x]
  
  reg_1 <- felm(trustlegisl~Dem+Rep+ownpres+ownpres:Rep+ownhouse+ownhouse:Rep+ownsenate+ownsenate:Rep++realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est_pres <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_levelpip_est_senate <- as.data.frame(summary(reg_1)$coefficients)["ownsenate",]
  fInfo_1_levelpip_est_house <- as.data.frame(summary(reg_1)$coefficients)["ownhouse",]
  fInfo_1_asymm_est_pres <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  fInfo_1_asymm_est_senate <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownsenate",]
  fInfo_1_asymm_est_house <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownhouse",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est_pres,fInfo_1_levelpip_est_senate,fInfo_1_levelpip_est_house,fInfo_1_asymm_est_pres,fInfo_1_asymm_est_senate,fInfo_1_asymm_est_house))
  fInfo_all$Effect <- c("Main effect President","Main effect Senate","Main effect House","Asymmetry President","Asymmetry Senate","Asymmetry House")
  fInfo_all$Government <- ifelse(grepl("Senate",fInfo_all$Effect),"Senate",
                                 ifelse(grepl("House",fInfo_all$Effect),"House","President"))
  fInfo_all$Effect <- gsub(" Senate| House| President","",fInfo_all$Effect)
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_decl_congress_conflegisl <- do.call(rbind.data.frame, lYearSens_decl_congress_conflegisl)

colnames(fInfo_all_decl_congress_conflegisl) <- c("Estimate","SE","T_val","P_val","Effect","Government", "MinYear")

fInfo_all_decl_congress_conflegisl <- as.data.table(fInfo_all_decl_congress_conflegisl)[,CI_lower := Estimate-1.96*SE]
fInfo_all_decl_congress_conflegisl <- as.data.table(fInfo_all_decl_congress_conflegisl)[,CI_upper := Estimate+1.96*SE]
fInfo_all_decl_congress_conflegisl <- fInfo_all_decl_congress_conflegisl[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]



pdf("Graphs/Window_estimates_decl_congress_conflegisl.pdf")
ggplot(fInfo_all_decl_congress_conflegisl, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Government) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()




#Panel D
#alignment on conf in congress, shifting window
lYearSens_shift_congress_conflegisl <- lapply(vMinYear, function(x){
  fData_rest <- fData[year%in%x:(x+iWindowLength)]
  
  reg_1 <- felm(trustlegisl~Dem+Rep+ownpres+ownpres:Rep+ownhouse+ownhouse:Rep+ownsenate+ownsenate:Rep++realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
  fInfo_1_levelpip_est_pres <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
  fInfo_1_levelpip_est_senate <- as.data.frame(summary(reg_1)$coefficients)["ownsenate",]
  fInfo_1_levelpip_est_house <- as.data.frame(summary(reg_1)$coefficients)["ownhouse",]
  fInfo_1_asymm_est_pres <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
  fInfo_1_asymm_est_senate <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownsenate",]
  fInfo_1_asymm_est_house <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownhouse",]
  
  fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est_pres,fInfo_1_levelpip_est_senate,fInfo_1_levelpip_est_house,fInfo_1_asymm_est_pres,fInfo_1_asymm_est_senate,fInfo_1_asymm_est_house))
  fInfo_all$Effect <- c("Main effect President","Main effect Senate","Main effect House","Asymmetry President","Asymmetry Senate","Asymmetry House")
  fInfo_all$Government <- ifelse(grepl("Senate",fInfo_all$Effect),"Senate",
                                 ifelse(grepl("House",fInfo_all$Effect),"House","President"))
  fInfo_all$Effect <- gsub(" Senate| House| President","",fInfo_all$Effect)
  fInfo_all$MinYear <- x
  return(fInfo_all)
  
})

fInfo_all_shift_congress_conflegisl <- do.call(rbind.data.frame, lYearSens_shift_congress_conflegisl)

colnames(fInfo_all_shift_congress_conflegisl) <- c("Estimate","SE","T_val","P_val","Effect","Government", "MinYear")

fInfo_all_shift_congress_conflegisl <- as.data.table(fInfo_all_shift_congress_conflegisl)[,CI_lower := Estimate-1.96*SE]
fInfo_all_shift_congress_conflegisl <- as.data.table(fInfo_all_shift_congress_conflegisl)[,CI_upper := Estimate+1.96*SE]
fInfo_all_shift_congress_conflegisl <- fInfo_all_shift_congress_conflegisl[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]



pdf("Graphs/Window_estimates_shift_congress_conflegisl_conflegisl.pdf")
ggplot(fInfo_all_shift_congress_conflegisl, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
  geom_point() +
  geom_smooth(se=F) +
  geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
  geom_hline(yintercept = 0)+
  facet_wrap(~Government) +
  theme_bw() +
  labs(fill="",col="")+
  xlab("First year in sample")+
  theme(text=element_text(size=18),legend.position = "bottom")
dev.off()










############# Figure A10: Evolution of president-in-power effect, shifting window of different lengths ################

#10 years problematic because there is a 10year period without a turnover
#estimate rolling window for each option between 10 and 20.
lapply(10:20, function(iWindowLength_diff){
  lYearSens_shift_diffyears <- lapply(vMinYear, function(x){
    fData_rest <- fData[year%in%x:(x+iWindowLength_diff)]
    
    reg_1 <- felm(trustfed~Dem+Rep+ownpres+ownpres:Rep+realinc|year+age+ female+ wrkstat+ degree+ race +relig,weights = fData_rest$wtssall, data=fData_rest)
    fInfo_1_levelpip_est <- as.data.frame(summary(reg_1)$coefficients)["ownpres",]
    fInfo_1_asymm_est <- as.data.frame(summary(reg_1)$coefficients)["Rep:ownpres",]
    
    fInfo_all <- as.data.frame(rbind(fInfo_1_levelpip_est,fInfo_1_asymm_est))
    fInfo_all$Effect <- c("Main effect","Asymmetry")
    fInfo_all$MinYear <- x
    return(fInfo_all)
    
  })
  
  fInfo_all_shift_diffyears <- do.call(rbind.data.frame, lYearSens_shift_diffyears)
  
  colnames(fInfo_all_shift_diffyears) <- c("Estimate","SE","T_val","P_val","Effect","MinYear")
  
  fInfo_all_shift_diffyears <- as.data.table(fInfo_all_shift_diffyears)[,CI_lower := Estimate-1.96*SE]
  fInfo_all_shift_diffyears <- as.data.table(fInfo_all_shift_diffyears)[,CI_upper := Estimate+1.96*SE]
  fInfo_all_shift_diffyears <- fInfo_all_shift_diffyears[,Effect := factor(Effect,levels=c("Main effect","Asymmetry"))]
  
  tryCatch({dev.off},error=function(e){})
  ggplot(fInfo_all_shift_diffyears, aes(x=MinYear,y=Estimate,col=Effect,fill=Effect)) +
    geom_point() +
    geom_smooth(se=F) +
    geom_ribbon(aes(ymin=CI_lower,ymax=CI_upper),alpha=0.1)+
    geom_hline(yintercept = 0)+
    # facet_wrap(~Effect) +
    theme_bw() +
    labs(fill="",col="")+
    xlab("First year in sample")+
    theme(text=element_text(size=18),legend.position = "bottom")
  ggsave(paste0("Graphs/Window_estimates_shift_diffyears_",iWindowLength_diff,".pdf"))
  
})

############ In-text: Negative weights and alternative diff-in-diff  ##############

#examine whether there are negative weights (there are not)
bacon_weights <- twowayfeweights(df=fData, 
                                 Y = "trustfed",
                                 G = "Party_cat",
                                 T = "year",
                                 D = "ownpres" ,
                                 cmd_type = "feTR")


# Create a tidier for "multiplegt" objects
tidy.did_multiplegt = function(x, level = 0.95) {
  ests = x[grepl("^placebo_|^effect|^dynamic_", names(x))]
  ret = data.frame(
    term      = names(ests),
    estimate  = as.numeric(ests),
    std.error = as.numeric(x[grepl("^se_placebo|^se_effect|^se_dynamic", names(x))]),
    N         = as.numeric(x[grepl("^N_placebo|^N_effect|^N_dynamic", names(x))])
  ) |>
    # For CIs we'll assume standard normal distribution
    within({
      conf.low  = estimate - std.error*(qnorm(1-(1-level)/2))
      conf.high = estimate + std.error*(qnorm(1-(1-level)/2))
    })
  return(ret)
}

#use did_multiplegt to compare to main effect (similar effect size)
reg_multiplegt <- did_multiplegt(df = fData,
                                 Y = "trustfed",
                                 G = "Party_cat",
                                 T = "year",
                                 D = "ownpres",
                                 brep      = 20                  # no. of bootstraps (required for SEs)
)

tidy.did_multiplegt(reg_multiplegt)









